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ABSTRACT 



We consider the nonhnear outcome of gravitational instability in optically- 
thick disks with a realistic cooling function. We use a numerical model that is 
local, razor-thin, and unmagnetized. External illumination is ignored. Cooling is 
calculated from a one-zone model using analytic fits to low temperature Rosseland 
mean opacities. The model has two parameters: the initial surface density So 
and the rotation frequency VL. We survey the parameter space and find: (1) The 
disk fragments when {{tc))Q ~ 1, where ((tc)) is an effective cooling time defined 
as the average internal energy of the model divided by the average cooling rate. 
This is consistent with earlier results that used a simplified cooling function. (2) 
The initial cooling time Tco for a uniform disk with Q = 1 can differ by orders of 
magnitude from ((tc)) in the nonlinear outcome. The difference is caused by sharp 
variations in the opacity with temperature. The condition t^o^ ~ 1 therefore 
does not necessarily indicate where fragmentation will occur. (3) The largest 
difference between {{tc)) and Tco is near the opacity gap, where dust is absent 
and hydrogen is largely molecular. (4) In the limit of strong illumination the 
disk is isothermal; we find that an isothermal version of our model fragments for 
Q < lA. Finally, we discuss some physical processes not included in our model, 
and find that most are likely to make disks more susceptible to fragmentation. 
We conclude that disks with {{Tc))il < 1 do not exist. 

Subject headings: accretion, accretion disks, solar system: formation, galaxies: 
nuclei 



1. Introduction 



The outer regions of accretion disks in both active galactic nuclei (AGN) and young stel- 
lar objects (YSO) are close to gravitational instability (for a review see, for AGN: Shlosman 
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et al. 1990; YSOs: Adams & Lin 1993). Gravitational instability can be of central impor- 
tance in disk evolution. In some disks, it leads to the efficient redistribution of mass and 
angular momentum (e.g. Larson 1984; Laughlin & Rozyczka 1996; Gammie 2001). In other 
disks, gravitational instability leads to fragmentation and the formation of bound objects. 
This may cause the truncation of circumnuclear disks (Goodman 2003), or the formation of 
planets (e.g. Boss 1997, and references therein). 

We will restrict attention to disks whose potential is dominated by the central object, 
and whose rotation curve is therefore approximately Keplerian. Gravitational instability to 
axisymmetric perturbations sets in when the sound speed Cg, the rotation frequency Vl, and 
the surface density S satisfy 

Q=^< QcrU ^ 1 (1) 

(Toomre 1964; Goldrcich & Lynden-Bell 1965). Here Qcrit = 1 ^oi a "razor-thin" (two- 
dimensional) fluid disk model of the sort we will consider below, and Qcrit = 0.676 for a 
finite-thickness isothermal disk (Goldrcich & Lynden-Bell 1965). ^ The instability condition 
(1) can be rewritten, for a disk with scale height H ~ Cg/fl, around a central object of mass 

Mdisk > -M„ (2) 
r 

where M^isk = vrr^E. For YSO disks H/r ~ 0.1 and thus a massive disk is required for 
instability. AGN disks are expected to be much thinner. The instability condition can be 
rewritten in a third, useful form if we assume that the disk is in a steady state and its 
evolution is controlled by internal ("viscous") transport of angular momentum. Then the 
accretion rate M = Snac^H/Q, where a < 1 is the usual dimensionless viscosity of Shakura 
& Sunyaev (1973), and 

M>5|i=7,><10-„(^;M,yr' (3) 

implies gravitational instability (e.g. Shlosman et al. (1990)). Disks dominated by external 
torques (e.g. a magnetohydrodynamic [MHD] wind) can have higher accretion rates (but 
not arbitrarily higher; see Goodman 2003) while avoiding gravitational instability. 

For a young, solar-mass star accreting from a disk with a = 10~^ at lO~^M0yr~^, 
equation (3) implies that instability occurs where the temperature drops below 17 K. Disks 
may not be this cold if the star is located in a warm molecular cloud where the ambient 



^For global models with radial structure, nonaxisymmetric instabilities typically set in for slightly larger 
values of Q (see Boss 1998 and references therein). 
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temperature is greater than 17 K, or if the disk is bathed in scattered infrared hght from 
the central star (although there is some evidence for such low temperatures in the solar 
nebula, e.g. Owen ct al. 1999). If the vertically-averaged value of a is small and internal 
dissipation is confined to surface layers, as in the layered accretion model of Gammie (1996), 
then instability can occur at higher temperatures, although equation (2) still requires that 
the disk be massive. 

AGN disk heating is typically dominated by illumination from a central source. The 
temperature then depends on the shape of the disk. If the disk is flat or shadowed, however, 
and transport is dominated by internal torques, one can apply equation (3). For example, 
in the nucleus of NGC 4258 (Miyoshi et al. 1995) the accretion rate may be as large as 
lO^^Mgyr"^ (Lasota et al. 1996; Gammie et al. 1999). Equation (3) then implies that 
instability sets in where T < 10^(q;/10~^) K. If the disk is illumination-dominated then Q 
fluctuates with the luminosity of the central source. 

In a previous paper (Gammie 2001), one of us investigated the effect of gravitational 
instability in cooling, gaseous disks in a local model. A simplifled cooling function A was 
employed in these simulations, with a flxed cooling time Tco'- 

A=-f. (4) 

' CO 

where V = the internal energy per unit area. Disk fragmentation was observed for Tco^ ^ 3. 
The purpose of this paper is to investigate gravitational instability in a local model with 
more realistic cooling. 

Several recent numerical experiments have included cooling, as opposed to isothermal 
or adiabatic evolution, and we can ask whether these results are consistent with Gam- 
mie (2001). Nelson et al. (2000) studied a global two-dimensional (thin) SPH model in 
which the vertical density and temperature structure is calculated self-consistently and each 
particle radiates as a blackbody at the surface of the disk. The initial conditions at a ra- 
dius corresponding to the minimum initial value of Q (~ 1.5) for these simulations were 
Eo ~ 50gcm~^,Q Ri S X 10~^°s~^; the initial cooling time under these circumstances is 
Tco ~ 250 so fragmentation is not expected and is not observed. 

Durisen et al. (2001) consider a global three dimensional (3D) Eulerian hydrodynamics 
model in which the volumetric cooling rate varies with height above the midplane so as to 
preserve an isentropic vertical structure. The coohng time is flxed at each radius. Their 

cooling time > 1011"^ at all radii, so fragmentation is not expected based on the criterion of 
Gammie (2001). The simulations show structure formation due to gravitational instabilities 
but not fragmentation. 
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Rice et al. (2003) consider a global 3D SPH model with a cooling time that is a fixed 
multiple of n~^(r). They find that their disk fragments when Tco ~ 3Q~^ and M^isk = 
O.IM*. For a more massive disk {M^isk = 0.25M^,), fragmentation occurred at somewhat 
higher cooling times [t^o ~ 10^2^^). This is effectively a global generalization of the local 
model problem considered by Gammie (2001). The fact that the results are so consistent 
suggests that the local, thin approximation used in Gammie (2001) and here give a reasonable 
approximation to a global outcome. 

Mayer et al. (2002) consider a global three dimensional SPH model of a circumstellar 
disk. Explicit cooling is not included, but the equation of state switches from isothermal 
to adiabatic when gravitational instability begins to set in. This is designed to account for 
the inefficient cooling of dense, optically thick regions. Fragmentation is observed. Realistic 
cooling can have a complex influence on disk evolution, and it is not clear that switching 
between isothermal and adiabatic behavior "brackets" the outcomes that might be obtained 
when full cooling is used. 

Other notable recent work, such as that by Boss (2002), includes strong radiative heating 
in the sense that the effective temperature of the external radiation fleld T^^r is comparable 
to or larger than the disk midplane temperature Tc- In the limit that Tirr <S Tc we recover 
the limit considered here and in Gammie (2001); in the limit that Ti^r ^ the disk is 
effectively isothermal. 

The plan of this paper is as follows. In §2 we describe the model, with a detailed 
description of the cooling function given in §3. The results of numerical experiments are 
described in §4. Conclusions are given in §5. 

2. Model 

The model we use here is identical to that used in Gammie (2001) in every respect 
except that we use a more complicated cooling function. To make the description more 
self-contained, we summarize the basic equations of the model here. The model is local, in 
the sense that it considers a region of size L where L/tq and To is a fiducial radius. We 

use a local Cartesian coordinate system x = r — and y = [(f) — Qt)ro, where r, are the 
usual cylindrical coordinates and Q is the orbital frequency at Tq. The model is also thin in 
the sense that matter is confined entirely to the plane of the disk. 

Using the local approximation one can perform a formal expansion of the equations of 
motion in the small parameter L/vo- The resulting equations of motion read, where v is the 
velocity, P is the (two-dimensional) pressure, and 4> is the gravitational potential with the 
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time-steady axisymmetric component removed: 



^ = - 2ne^ X V + Sifxe^ - V0. (5) 

For constant pressure and surface density, v = —^flxey is an equilibrium solution to the 
equations of motion. This linear shear flow is the manifestation of differential rotation in 
the local model. 

The equation of state is 

p = (7 - m (6) 

where P is the two-dimensional pressure and U the two-dimensional internal energy. The 
two-dimensional (2D) adiabatic index 7 can be mapped to a 3D adiabatic index F in the low- 
frequency (static) limit. For a non-self-gravitating disk 7 = (3F — 1)/(F + 1) (e.g. Goldreich, 
Goodman, & Narayan 1986; Ostriker, Shu, & Adams 1992). For a strongly self-gravitating 
disk, one can show that 7 = 3 — 2/F. We adopt T — 7/5 throughout, which yields 7 = 11/7. 

The internal energy equation is 

dU 

_+V-(C/v) = -PV-v-A, (7) 

where A = A(S, U, fl) is the cooling function, fully described below. Notice that there is no 
heating term; heating is due solely to shocks. Numerically, entropy is increased by artificial 
viscosity in shocks. 

The gravitational potential is determined by the razor-thin disk Poisson equation: 

V'^(f) = A7rGE6{z). (8) 
For a single Fourier component of the surface density Ek this has the solution 

= _^Eke*-^-l'=^l. (9) 



A finite-thickness disk has weaker self-gravity, but this does not qualitatively change the 
dynamics of the disk in hnear theory (Goldreich & Lynden-Bell 1965). 

We integrate the governing equations using a self-gravitating hydrodynamics code based 
on ZEUS (Stone & Norman 1992). ZEUS is a time-explicit, operator-split, finite-difference 
method on a staggered mesh. It uses an artificial viscosity to capture shocks. Our imple- 
mentation has been tested on standard linear and nonlinear problems, such as sound waves 
and shock tubes. We use the "shearing box" boundary conditions, described in detail by 
Hawley et al. (1995), and solve the Poisson equation using the Fourier transform method. 
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modified for the shearing box boundary conditions. See Gammie (2001) for further details 
on boundary conditions, numerical methods and tests. 

The numerical model is always integrated in a region of size L x L at a numerical 
resolution oi N x N. In linear theory the disk is most responsive at the critical wavelength 
2c3/GSo.^ We have checked the dependence of the outcome on L and found that as long as 

L > 2(?g/GTio the outcome does not depend on L. We have also checked the dependence of 
the outcome on and found that the outcome is insensitive to A^, at least for the models 
with N > 256 that we use. 



3. Cooling Function 

Our cooling function is determined from a one-zone model for the vertical structure of 
the disk. The disk cools at a rate per unit area 

A = 2aT^, (10) 

which defines the effective temperature Tg. The cooling function depends on the heat content 
of the disk and how that content is transported from the disk interior to the surface: by 
radiation, convection, or perhaps some more exotic form of turbulent transport such as MHD 
waves. Low temperature disks are expected to be convectively unstable (e.g. Cameron 1978; 
Lin & Papaloizou 1980). Cassen (1993) has argued, however, that the radiative heat flux 
in an adiabatically-stratified disk is comparable to the heat dissipated by turbulence (in 
an a-disk model), suggesting that convection is incapable of radically altering the vertical 
structure of the disk. We will consider only radiative transport. 

If the disk is optically thick in the Rosseland mean sense, so that radiative transport 
can be treated in the diffusion approximation, then (Hubeny 1990) 

8 T'^ 

Tt = (11) 

where r is the Rosseland mean optical depth and is the central temperature. We will 
assume that Tc ~ T, where 

T=^^^, (12) 

and 

c^ = 7(7-l)|, (13) 



^The wavelength corresponding to the minimum in the dispersion relation for axisymmetric waves. 
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which follows from the equation of state and the assumption that the radiation pressure 
is small (we have verified that this is never seriously violated). Here ks is Boltzmann's 

constant, rrip is the proton mass, and n is the mean mass per particle, which we have set to 
2.4 in models with initial temperature below the boundary between the grain-evaporation 
opacity and molecular opacity and p = 0.6 in models with initial temperature above the 
boundary. 

The optical depth is 

poo 

T= dzK{p„T,)p, (14) 

^0 

where k is the Rosseland mean opacity, pz and are local density and temperature, and z 
is the height above the midplane. Following the usual one-zone approximation, 

POO 

/ dzK{pz,Tz)pz ^ Hk{p,T)p (15) 
Jo 

where the overbar indicates a suitable average and H ^ Cs{T)/Q is the disk scale height (we 
ignore the effects of self-gravity on the disk scale height, which is vahd when locally Q > 1). 
Taking T and p ^ T,/{2H) then gives a final, closed expression for A. 

We have adopted the analytic approximations to the opacities provided by Bell & Lin 
(1994). These opacities are dominated by, in order of increasing temperature: grains with ice 
mantles, grains without ice mantles, molecules, H~ scattering, bound- free/free-free absorp- 
tion and electron scattering. The molecular opacity regime is commonly called the opacity 
gap; it is too hot for dust, but too cold for H~ scattering to contribute much opacity. The 
opacity can be as much as 4 orders of magnitude smaller than the ~ 5gcm^^ typical of the 
dust-dominated opacity regime. It turns out that this feature plays a significant role in the 
evolution of gravitationally-unstable disks. 

To sum up, the cooling function is 

A(E,C/,(]) = i^^. (16) 
o r 

For a power- law opacity of the form k — kqp"'T'^, this implies that 

A J2-5-3a/2+bjj4+a/2-b _ ^^7^ 

From this it follows that the cooling time Tc = U/A scales as 

Tc ~ J]5+3a/2-6fj-3-a/2+&_ ^j^g^ 

If the disk evolves quasi-adiabatically (as it does if the cooling time is long compared to the 
dynamical time) then t/ ~ E'*^ and 

Tc ~ (19) 
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Table 1 gives a list of values for this scaling exponent for our nominal value of 7 = 11/7. 
Notice that, when ice grains or metal grains are evaporating, and in the bound- free/free- free 
opacity regime, cooling time decreases as surface density increases. 

Our cooling function is valid in the limit of large optical depth (r » 1). Since the disk 
becomes optically thin at some locations in the course of a typical run, we must modify this 
result so that the cooling rate does not diverge at small optical depth. A modification that 
produces the correct asymptotic behavior is 



A ^ '^oT-—, (20) 



This interpolates smoothly between the optically-thick and optically-thin regimes and is 
proportional to the (Rosseland mean) optical depth in the optically-thin limit. While it 
would be more physically sensible to use a Planck mean opacity in the optically-thin limit, 
usually the optically-thin regions contain little mass so their cooling is not energetically 
significant. An exception is in the opacity gap, where even high density regions become 
optically thin. 

Our simulations begin with S and U constant. The velocity field is perturbed from 
the equilibrium solution to initiate the gravitational instability. The initial velocities are 
V.J. = 6v.j., Vy = —^Qx + 6vy, where 5v is a Gaussian random field of amplitude {Sv'^)/c'j, = 0.1. 
The power spectrum of perturbations is white noise {vl ~ in a band in wavcnumber 
kcrit/4: < \k\ < 4:kcrit surrounding the minimum kcru = ^/{''^Q^) (with G = So = = 
1) in the density- wave dispersion relation. We have checked in particular cases that for 
10~^ < {5v'^)/c1 < 10 the outcome is qualitatively unchanged. This is expected because 
disk perturbations (unlike cosmological perturbations) grow exponentially and the initial 
conditions are soon forgotten. 

Excluding the initial velocity field, the initial conditions for a spatially-uniform disk 
consist of three parameters: Eq, Ug, and Q. We fix Q = 1, leaving two degrees of freedom. In 
models with simple, scale- free coohng functions such as that considered by Gammie (2001), 
these degrees of freedom remain and can be scaled away by setting G = T^^ = Q = 1. That 
is, there is a two-dimensional continuum of disks (with varying values of Eq and fl, but the 
same value of Q) that are described by a single numerical model. 

The opacity contains definite physical scales in density and temperature. The realistic 
cooling function considered here therefore removes our freedom to rescale the disk surface 
density and rotation frequency. That is, there is now a one-to-one correspondence between 
disks with fixed E and fl and our numerical models. 



The choice of E,, and Q as labels for the parameter space is not unique. Internally in the 
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code we fix the initial volume density (in gcm"^) and the initial temperature (in Kelvins). 
These choices are difficult to interpret, however, since they are tied to quantities that change 
over the course of the simulation; fl and the mean value of E do not. 

The coohng is integrated explicitly using a first-order scheme. The timestep is modified 
to satisfy the Courant condition and to be less than a fixed fraction of the shortest cooling 
time on the grid. We have varied this fraction and shown that the results are insensitive to 
it, provided that it is sufficiently small. 



4. Nonlinear Outcome 

4.1. Standard Run 

Consider the evolution of a single "standard" run, with Eq = 1.4 x lO^gcm"^ and 
n = l.lx lO^^sec-^ This corresponds to To = 1200 and r^o = 9.0 x 10^n-\ The model 
size is L = 320G'E„/f2^ and numerical resolution 1024^. The model initially lies at the lower 
edge of the opacity gap. 

The evolution of the kinetic, gravitational and thermal energy per unit area {{Ek), 
(Eg) and (Eth) respectively) normalized to G^E^/fi^,^ are shown in Figure 1. After the 
initial phase of gravitational instability the model settles into a statistically-steady, gravito- 
turbulent state. It does not fragment. Coohng is balanced by shock heating. Energy for 
driving the shocks is extracted from the shear fiow, and the mean shear fiow is enforced by 
the boundary conditions. 

The turbulent state transports angular momentum outward via hydrodynamic and grav- 
itational shear stresses. The dimensionless gravitational shear stress is 



1 



oo 



'''''''''' il^cj) J _ J' to ^^^^ 

where g is the gravitational acceleration, and the dimensionless hydrodynamic shear stress 
is 

(^hyd = 73TT3f (22) 



where () denote a spatial average. Figure 2 shows the evolution of {ag^av) and {ahyd) in the 
standard run. Averaged over the last 230Q~^ of the run, {{ahyd)) — 0.0079, {{agrav)) ~ 0.017, 



^The natural unit that can be formed from G, S and CI. 
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and so the total dimensionless shear stress is {{a)) — 0.025, where (()) denote a space and 
time average. 

The mean stabihty parameter (Q) = {cs)0./ttG{'E) averages 1.86 over the last 230^1"^ 
of the run. Because the temperature and surface density vary strongly, other methods of 
averaging Q will give different results. 

Figure 3 shows a snapshot of the surface density att — 50Q~^. The structure is similar to 
that observed in Gammie (2001), with trailing density structures. The density structures are 
stretched into a trailing configuration by the prevailing shear flow. Their scale is determined 
by the disk temperature and surface density rather than the size of the box (see Gammie 
2001). 

4.2. Varying So and 

We now turn to exploring the two-dimensional parameter space of models. First consider 
a series of models with the same initial central temperature, but with varying Tco- As Tco 
is lowered the time-averaged gravitational potential energy per unit area {{Eg)) increases 
monotonically in magnitude. The gravito-turbulent state becomes more extreme, with larger 
((a)), larger perturbed velocities, and larger density contrasts. Eventually a threshold is 
crossed and the disk fragments. 

Fragmentation is illustrated in Figure 4, which shows a snapshot from a run with Eo = 
6.6 X lO^gcm-^ ^] = 5.4 X IQ-^sec-^ This corresponds to = 1200, t^^ = 0.025n-\ 
The run has numerical resolution 256^ and L = SOGSo/fi^. The largest bound object in the 
center of the figure was formed from the collision and coalescence of several smaller bound 
objects. A snapshot of the optical depth at the same point in the simulation is given in 
Figure 5. For each snapshot, red indicates high values of the mapped variable and blue 
indicates low values. Much of the disk is optically thick, but most of the low density regions 
are optically thin in the Rosseland mean sense. 

Lowering Tco sufficiently always leads to fragmentation. We have surveyed the parameter 
space of Q, and Eg to determine where the disk begins to fragment. Each model was run 
to 100Q~^^ Figures 6 and 7 summarize the results. Two heavy sohd lines are shown on 
each diagram. The upper line shows the most rapidly cooling simulations that show no signs 
of gravitational fragmentation {nonfragmenfMion point). Quantitatively, we define this as 
the point at which the time-averaged gravitational potential energy per unit area is equal 



In four cases we had to run the simulation longer to get converged results. 
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to — SG^Ep/Q^.^ The lower line shows the most slowly cooling simulations to show definite 
fragmentation {fragmentation point) . Quantitatively, we define this as the point at which the 
gravitational potential energy per unit area is equal to — 300G^S^/r2^ at some point during 
the run.^ Figure 6 shows the data in the po, Tq plane, while Figure 7 shows the results in the 
So, Vt plane. The light contours are lines of constant Tco- 

The transition from persistent, gravito-turbulcnt outcomes to fragmentation is gradual 
and statistical in nature. Figure 8 shows the gravitational potential energy per unit area 
in the transition region for a series of runs with = 1200 K. The abscissa is labeled with 
the initial cooling time Tcofi. There is a gradual, approximately logarithmic increase in the 
magnitude of {{Eg)) as Tco decreases. Runs in this region exhibit the transient formation 
of small bound objects, which might collapse if additional physics (e.g. the effects of MHD 
turbulence) were included in the model. Eventually —{{Eg)) begins to increase dramatically, 
and we define the transition point as the beginning of this steep increase in gravitational 
binding energy. 

Figure 9 shows the run of Tco^ for the fragmentation point, transition point, and non- 
fragmentation point as a function of Tg. It is surprising that a disk can begin to exhibit 
signs of gravitational collapse for Tco^ as large as 10^, and evade collapse for Tco^ as small as 
0.02. A naive application of the results of Gammic (2001) would suggest that fragmentation 
should occur for Tco^ ^ 3. Evidently this estimate can be off by orders of magnitude, with 
the largest error for Tg fti 10^ K, just below the opacity gap. 

The physical argument for fragmentation at short coohng times is as follows (e.g. Shlos- 
man et al. 1990). Thermal energy is supplied to the disk via shocks. Strong shocks occur 
when dense clumps collide with one another; this occurs on a dynamical timescale ~ 
If the disk cools itself more rapidly then shock heating cannot match cooling and fragmen- 
tation results. This argument is apparently contradicted by Figure 9. The resolution lies 
in finding an appropriate definition of cooling time. The disk loses thermal energy on the 
effective cooling timescale 

(W>- - (23) 

Figure 10 shows the run of ((tc)) at the fragmentation, transition, and non- fragmentation 
points. Evidently ((tc)) at transition lies between and 10^2"^. Figure 11 shows the run 
of Tco s-nd ((tc)) on the transition line. Just below the opacity gap they differ by as much as 



^—3 is the potential energy per unit area of a wave at the eritieal wavelength in a Q = 1 disk with 
— VS/tt. No bound objects are observed throughout the duration of these runs. 

^These runs exhibit bound objects that persist for the duration of the run. 
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four orders of magnitude. 

Why do Tco and ((tc)) differ by such a large factor? The answer is related to the existence 
of sharp variations in opacity with temperature. Consider a disk near the lower edge of the 
opacity gap. Once gravitational instability sets in, fluctuations in temperature move parts 
of the disk into the opacity gap. There, the opacity is reduced by orders of magnitude. Since 
the cooling rate for an optically thick disk is proportional to k~^, the cooling time drops by a 
similar factor. Relatively small variations in temperature can thus produce large variations 
in cooling rate. 

As in Gammie (2001), the result ((tc))Q > 1 also imphes a constraint on {{a)). Energy 
conservation implies that 

Inm^y)) = ((A)), (24) 

where W^y is the total shear stress (hydrodynamic plus gravitational). Equivalently, stress 
by rate-of-strain is equal to the dissipation rate. Using the definition of ((tc)), this implies 

((«)) = (7(7 -l)^^^((rc))) (25) 

Hence {{Tc))ft > 1 implies ((a)) ^ 1. Figure 12 shows {{a)) vs ((tc)) for a large number of 
runs plotted against equation (25). For small values of ((tc)) the numerical values lie below 
the line. These models are not in equilibrium (i.e., not in a statistically-steady gravito- 
turbulent state), so the time average used in equation (24) is not well defined. For larger 
values of ((tc)) numerical results typically (there is noise in the measurement of both ((a)) 
and ((tc)) because the time average is taken over a finite time interval) lie slightly above the 
analytic result. 

The bias toward points lying slightly above the line reflects the fact that {{a)) measures 
the rate of energy extraction from the shear while ((tc)) measures the rate at which that 
energy is transformed into thermal energy. If energy is lost, perhaps to numerical averaging 
at the grid scale, then more energy must be extracted from the shear flow to make up the 
difference. Overall, however, the agreement with the analytic result is good and demonstrates 
good energy conservation in the code. 

The relationship between ((tc)) and {{a)) is interesting but not particularly useful be- 
cause ((tc)) is no more readily calculated than ((a)); it depends on a complicated moment 
of the surface density and temperature. Only for constant cooling time have we been able 
to evaluate this moment analytically. 
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4.3. Isothermal Disks 

We have assumed that external illumination of the disk is negligible. This approximation 
is valid when the effective temperature Tj„. of the external irradiation is small compared to 
the central temperature of the disk. In the opposite limit, illumination controls the energetics 
of the disk and it is isothermal (if it is illuminated directly so that shadowing effects, such 
as those considered by Jang-Condell & Sasselov (2003) are negligible). 

It is therefore worth studying the outcome of gravitational instability in an isothermal 
disk. The isothermal disk model has a single parameter: the initial value of Q. We ran 
models with varying values of Q and with {5v'^)/cl — 0.1. We find that models with Q < 1.4 
fragment. 

It is likely that the mass of the fragments, etc., depends on how an isothermal disk 
becomes unstable. Rapid fluctuation of the external radiation fleld is likely to produce a 
different outcome than dimming on a timescale long compared to the dynamical time. 

5. Discussion 

Using numerical experiments, we have identified those disks that are likely to fragment 
absent external heating. Disks with effective cooling times ((tc)) < are susceptible to 
fragmentation. This is what one might expect based on the simple argument of Shlosman et 
al. (1990): if the disk cools more quickly than the self-gravitating condensations can coUide 
with one another, then those collisions (which occur on a timescale ~ Q~^) cannot reheat 
the disk and fragmentation is inevitable. But our results are at the same time surprising. 

The effective cooling time depends on the nonlinear outcome of gravitational instability. 
It depends on the cooling function, which in turn depends sensitively on E and U. Since E 
and U vary strongly over the disk once gravitational instability has set in, it is difficult to 
estimate ((tc)) directly. One might be tempted to estimate ((rc))(E, fl) ~ TcoiX'o, = 1), 
but our experiments show that this estimate can be off by as much as four orders of magni- 
tude. The effect is particularly pronounced near sharp features in the opacity. For example, 
consider a model initially located just below the opacity gap with Tco^ ^ 1- Gravitational 
instability creates dense regions with higher temperatures, where dust is destroyed. The 
result is rather like having to shed one's blanket on a cold winter morning: the disk loses its 
thermal energy suddenly. Pressure support is lost and gravitational collapse ensues. 

The difference between ((Tc)) and Tco{Q = 1) implies that a much larger region of the disk 
is susceptible to fragmentation than naive estimates based on the approximation ((tc)) ~ Tco 
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might suggest. For example, consider an equilibrium disk model with Q ^ 1 at small r. 
As r increases, Q declines. Eventually Q ~ 1 and gravitational instability sets in. There 
is then a range of radii where Q ~ 1, ((rc))fi > 1 and recurrent gravitational instability 
can transport angular momentum and prevent collapse. Generally speaking, however, the 
cooling time decreases with increasing radius. Eventually {{tc))Q ~ 1 and fragmentation 
cannot be avoided. By lowering our estimate of ((tc)), we narrow the range of radii over 
which recurrent gravitational instability can occur. 

The general sense of our result is that it is extremely difficult to prevent a marginally- 
stable, Q ~ 1, optically-thick disk from fragmenting and forming planets (in circumstellar 
disks) or stars (in circumstellar and circumnuclear disks). This is particularly true for disks 
with T ~ 10^ K, whose opacity is dominated by dust grains, i.e. disks whose temperature 
lies within a factor of several of the opacity gap. 

Our numerical model uses a number of approximations. First, our treatment is razor- 
thin, i.e. all the matter is in a thin slice ed, z — 0. The effect of finite thickness on linear 
stabihty has been understood since Goldreich & Lynden-Bell (1965): it is stabilizing because 
gravitational attraction of neighboring columns of disk is diluted by finite thickness. The 
size of the effect may be judged by the fact that Q — 0.676 is required for marginal stability 
of a finite-thickness, isothermal disk. 

The behavior of a finite-thickness disk in the nonlinear regime is more difficult to predict. 
Shocks will evidently deposit some of their energy away from the midplane, where it can 
be radiated away more quickly (because the energy is deposited at smaller optical depth - 
see Pickett et al. 2000). Radiative diffusion parallel to the disk plane (not included here) 
may enhance cooling of dense, hot regions. Both these effects arc destabilizing. Ultimately, 
however, a numerical study is required. This is numerically expensive: one must resolve the 
disk vertically, on the scale height H, and horizontally, at the critical wavelength 2tiQH . 

Second, we have ignored magnetic fields. While there may be astrophysical situations 
where cool disks have such low ionization that they are unmagnetized, most disks are likely 
to contain dynamically important magnetic fields that give rise to a dimensionless shear 
stress ((a)) > 0.01 (e.g. Hawley et al. 1995). These fields are likely to remove spin angu- 
lar momentum from partially collapsed objects, destabilizing them. Numerical experiments 
including both gravitational fields and magnetohydrodynamics are necessarily three dimen- 
sional (the instability of Balbus & Hawley (1991) requires 9^ 7^ 0), and are thus numerically 
expensive. 

Third, we have fixed 7 and ji for the duration of each simulation. This eliminates the soft 
spots in the equation of state associated with ionization of atomic hydrogen and dissociation 
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of molecular hydrogen. In these locations the three dimensional 7 dips below 4/3, which is 
destabilizing. 

Fourth, we have treated the physics of grain destruction and formation very simply. In 
using the Bell & Lin (1994) opacities we implicitly assume that grains reform in cooling gas 
on much less than a dynamical time. It is likely that grain re-formation will take some time 
(e.g. Hessman 1991) and this will further reduce the disk opacity and enhance fragmentation. 

Fifth, we have neglected the effects of illumination. In the limit of strong external illu- 
mination, i.e. when the effective temperature of the irradiation Tj^r- is large compared to the 
disk central temperature Tg, the disk is isothermal (here Tc is the temperature of a dense con- 
densation). We have carried out isothermal experiments and shown that, for initial velocity 
perturbations with {6v^)/c1 = 0.1, disks with Q < 1.4 fragment. Weaker illumination pro- 
duces a more complicated situation that we have not explored here. Illumination-dominated 
disks that become unstable presumably do so because the external illumination declines, and 
the rate at which the external illumination changes may govern the nonlinear outcome. 

We conclude that disks with {{tc))Q < 1 do not exist. Cooling in this case is so effective 
that fragmentation into condensed objects- stars, planets, or smaller accretion disks- is 
inevitable. 

As an example application of this result, consider the model for the nucleus of NGC 1068 
recently proposed by Lodato & Bertin (2003). Their model is an extended marginally-stable 
sclf-gravitating disk of the type investigated here and originally proposed by Goldrcich & 
Lyndcn-BcU (1965) for galactic disks and Paczynski (1978) for accretion disks, although their 
disk is sufficiently massive that it modifies the rotation curve as well. Based on their Figure 
3, at a typical radius of 0.5 pc, Eq ~ 10^ and Q ~ 10~^. According to our Figure 7 this disk 
is about 2 orders of magnitude too dense to avoid fragmentation. While it may be possible 
to avoid this conclusion by invoking strong external heating, the energy requirements are 
severe, as outlined in Goodman (2003). The disk proposed by Lodato & Bertin (2003) would 
therefore fragment into stars on a short timescale. 

This work was supported by NSF grant AST 00-03091, PHY 02-05155, and NASA grant 
NAG 5-9180. 
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Table 1. Scaling Exponent for Cooling Time as a Function of Surface Density 



Opacity Regime 


a 


b 


Exponent 


Ice grains 





2 


10/7 


Evaporation of ice grains 





-7 


-26/7 


Metal grains 





1/2 


4/7 


Evaporation of metal grains 


1 


-24 


-89/7 


Molecules 


2/3 


3 


52/21 


H~ scattering 


1/3 


10 


131/21 


Bound-free and free-free 


1 


-5/2 


-3/7 


Electron scattering 








2/7 



-19- 




Fig. 1. — Evolution of the kinetic, gravitational, and thermal energy per unit area, normal- 
ized to G^Ej^/r^^, in the standard run, which has L — 320GEo/n^, resolution 1024^, and 
Tco = 9.0 X lO^O-^ 
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Fig. 2. — Evolution of the gravitational (solid line) and hydrodynamic (dashed line) pieces 
of {a) in the standard run. 



Fig. 3. — Map of surface density at t = 50^2 ^ in the standard run. Blue is low density 
(0.2So) and yellow is high density (3So). 




Fig. 5. — Map of optical depth in a run with Tco = 0.025f2 ^. Black is low optical depth 
(10^^) and red is high optical depth (10^). 
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Fig. 6. — Location of the critical curves as a function of initial volume density and temper- 
ature (in cgs units). Each contour line is an order of magnitude change in Tco, solid/dotted 
lines indicating positive/negative integer values of log(Tco). 
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Fig. 7. — Location of the critical curves as a function of initial surface density and rota- 
tion frequency (in cgs units). Each contour line is an order of magnitude change in Tco, 
solid/dotted lines indicating positive/negative integer values of log (tco). The gap in the 
center of the plot is due to the discontinuous jump in the value of /i. 




Fig. 8. — Mean gravitational potential energy as a function of initial cooling time for a series 
of models with varying initial cooling time and Tq = 1200. 
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Fig. 9. — Initial cooling times at the points of non-fragmentation, fragmentation and tran- 
sition. 
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Fig. 10. — Effective cooling times at the points of non-fragmentation, fragmentation and 
transition. 
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Fig. 11. — Initial and effective cooling times at the transition between non-fragmentation 
and fragmentation. 
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Fig. 12. — Time- averaged shear stress vs. effective cooling time for a series of runs. The 
sohd hne shows the analytic result, based on energy conservation, from equation (25). 



